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Abstract 



In a previous work (Callegari and Yokoyama 2007, Cclest. Mech. Dyn. Astr. vol. 98), the main 
features of the motion of the pair Enceladus-Dione were analyzed in the frozen regime, i.e., without 
considering the tidal evolution. Here, the results of a great deal of numerical simulations of a pair 
of satellites similar to Enccladus and Dione crossing the 2:1 mean-motion resonance are shown. The 
resonance crossing is modeled with a linear tidal theory, considering a two-degrees-of-freedom model 
written in the framework of the general three-body planar problem. The main regimes of motion of the 
system during the passage through resonance are studied in detail. We discuss our results comparing 
them with classical scenarios of tidal evolution of the system. We show new scenarios of evolution of the 
Enceladus-Dione system through resonance not shown in previous approaches of the problem. 



1. Introduction 



Due to the proximity of the 2:1 mean-motion resonance, the pair Enceladus-Dione (hereafter denoted 
by E-D) is currently trapped in a libration state where the line of conjunctions between the satellites 
oscillates around the pericenter of Enccladus with period of 11.19 years and amplitude of 1.505° (Sinclair 
1972, 1983; Christou et al. 2007). The current resonant orbital configuration of the Saturnian satellites 
is not primordial, and probably was reached during the tidal-induced migration of the satellites due 
to forces acting on the planet and on the satellites (Goldreich 1965; Sinclair 1972, 1974; Sinclair 1983; 
Henrard and Lcmaitre 1983; Peale 1986, 1999, 2003; Meyer and Wisdom 2008). The formation of the 
current configuration of the system is not well understood yet in spite of much effort done in the works 
mentioned above. 

Several authors (e.g. Sinclair 1972; Henrard and Lemaitre 1983; Peale 1986) have studied the evolution 
of the E-D system through resonance with simplified versions of the general three-body problem given by 
one-degree-of-freedom models. The justification to use simplified models of resonance is based on the fact 
that the multiplet of resonances at 2:1 commensurability are well separated since Saturn's oblateness is 
large (see discussion in Peale 1986). Some works (Ferraz-Mello and Dvorak 1987; Callegari and Yokoyama 
2007 — hereafter denoted by CY2007; Meyer and Wisdom 2008), explore the dynamics of a pair of 
satellites similar to E-D with more general models of the three-body problem. In CY2007, a two-degrees- 
of-freedom model is developed to study in details the loci and the domain of the main and secondary 
resonances inside the 2:1 resonance. Only the conservative dynamics of the system in different portions 
of the phase space is studied in CY2007. Meyer and Wisdom (2008) investigate the past and future 
evolution of the pair Enceladus-Dione with a general model of the 2:1 resonance which also includes tidal 
dissipation on Saturn, Enceladus and Dione. 

Here, the non-conservative dynamics of Enccladus and Dione during the passage through 2:1 resonance 
is investigated with an extension of the previous model (CY2007), including now the time variations of 
the coefficients of the model due to tides on Saturn. We show the results of evolution through resonance 
considering initial parameters in a wide portion of the phase space. Our methodology is similar to that 
developed by Tittemore and Wisdom (1988, 1990) (hereafter denoted by TW1988, TW1989, TW1990, 
respectively) and Peale (1988, 1999). These authors investigated the dynamics of several resonances 
amongst the Uranian satellites. In particular, Peale (1988), TW1990 and Peale (1999) discussed the 
dynamics of the pair Ariel and Umbriel through the 2:1 resonance. TW1990 estimated the probabilities 
of capture and escape from the 2:1 resonance in a wide range of initial values of eccentricities. Here we 
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follow the methodology of TW1990 and apply it to the case of Enccladus and Dione. 

This paper has two main purposes: 1) to show a detailed description of the dynamics of a pair of 
regular satellites crossing the 2:1 mean- motion resonance in the planar case under the action of oblatcncss 
and perturbations of the central body; 2) to discuss the possibility to recover the past history of the pair 
E-D with our two-degrees-of-freedom model. 

We have divided the work as follows. The model and the initial parameters are given briefly in Section 
2. The results of a large deal of numerical simulations are shown in Section 3, where a classification of 
the types of the motion during the passage of the resonance is also given. In Section 4, we discuss the 
results obtained in Section 3 comparing them with those obtained with integrable models and most recent 
works. Section 5 is devoted to general conclusions. 



2. Modeling the passage through the 2:1 resonance 
2.1 Model 

In this work we use the averaged-planar model developed in CY2007. The model is based on the classical 
Laplacian expansions of disturbing function, and is constructed in the framework of the general three- 
body problem. The secular effects due to the perturbations of the non-spherical shape of Saturn arc 
included considering J 2 and J4 oblatcncss terms. Only the main resonant and secular terms are kept in 
the expression of the disturbing function. 

Denoting by mi, aj, n%, ej, Ai, vji the values of mass, semi- major axis, mean motion, eccentricity, 
mean longitude and the longitude of pericenter of the Enceladus (subscript E) and Dione (subscript D), 
respectively, we consider the following set of resonant variables 

Ie = Le — Ge, oe = 2 Ad — Xe — vje', 

Id = Ld — Gd, cfd = 2 Ad — Xe — zud; 

T e = L E + (Ie + Id), X e ; (1) 

T_D = Ld — 2(Ie + Id), ^d, 



where L { = fi iy fjl~a~ l , d = L^l - ef, fa = , m =G(M + m 4 ), M is the Saturn mass, 0=1293.4686 

is the gravitational constant in units of day, equatorial radius of Saturn (R e ) and M, whose values are 
listed in Table I. 

In terms of the non-singular variables Xi = cosci, yi = y/2li sin ai, and neglecting the constant 
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terms, the expanded Hamiltonian can be written as 

H = A(x%+y% + x 2 D +y 2 D ) 

+B(x 2 E + y% + x 2 D + y 2 D f 

+C{x% + y%) + D{x 2 D + y 2 D ) + E{x E x D + y E y D ) (2) 
+Fx E + Ix D 

+K(x E x D - y E yo) + R{x 2 E - y%) +S(x 2 D - y 2 D ), 

where A — F, I, K , R and S are the coefficients of the Hamiltonian which depend on the masses, J2, J4, 
eccentricities, semi-major axes, and the ratio a = ^ through combinations of the Laplace coefficients. 
The expressions and numerical values of the coefficients are given in Table II (first and third columns). 
They are calculated in the exact 2:1 commensurability (see Table I). In terms of resonant variables defined 
by Eq. ([T]), the coefficients are multiplied by the quantities given in the second column in Table II. Note 
that since the averaged Hamiltonian ^ is cyclic in and X E , T E and T^j are constants of motion, and 
the dynamical system has two degrees of freedom. 
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Table I: Model parameters and tidal constants. 



Sati 



Enceladus 



Dionc 



J2-, Ji a 
Mass a 
Radius (km) a - fc 
Q: dissipation function c 
k2'- dynamical Love number c 

Eccentricity 
Semi-major axes in the exact 
2:1 commensurability (R e ) 
Initial semi-major axis^ 
Current semi-major axis b 
Current eccentricity 6 ' 3 



^0 



0.0162907, -0.00093583 
3498.96" 1 Mr. 
60268 
>18,000 
0.34 



1.901 x 10~ 7 A/ 19.275 x 10" 7 M 



252.3 
100 d 
0.00075 d 


3.938806 

3.938579 
3.955812 
0.0045 



562.5 
100 e 
0.00075 6 


6.252469 

6.252288 
6.266127 
0.0022 



a Jacobson et al. (2006). M@ is the solar mass. 

b http://ssd.jpl.nasa.gov/horizons.html January 2007. 



c Murray and Dermott (1999). 

d In Meyer and Wisdom (2008), Q E = 20 and k 2E = 0.0018 are also adopted. 

e Supposing the same value as the Enceladus one. 

' Values considered in the calculations of the initial 5 parameter. 

9 See also Sinclair (1972), Ferraz-Mello (1985). 



4 



Table II: Expressions of the coefficients of the Hamiltonian and integrals I\ and their 
numerical values. 



Coefficient 



Terms in Hamiltonian Numerical Value 



2r E rj 



C=\k 12 [: 

D = \k 12 {^[- 



A = 
B = 

4a 

-4a| 



hGM 



\GM 



k E 
k E 

[r% ' 

2 da PpUD 

2 da PnHD 
'W-daW^ 



2k L 
4k L 



^E I 

D D 



D 



} 



2c "1" 2j4c 
7j2rf + Tj4d 



E = -- 



ck 12 



2^^ 
efci2 



_ hk 12 



R = 
S 



/fcl2 

sr B riL 



gfci: 
8T? 



2(/ B + Id) 




-5.2838 x 10~ 7 


4(Ie + Id) 2 




-147057.5727 


2I e 




-8.043543 x 10~ 3 


2I D 




-5.1166716 x 10~ 3 


2^I E I D cos(ct e - 


o\d) 


8.998116 x 10~ 7 


\ / 2I E COSCTe 




2.4498297 x 10~ 8 


\J2Id cos (Td 




-6.0934146 x 10~ 9 


2^I E I D cos(a E + 




7.7631926 x 10~ 6 


21 E cos 2(7 £ 




-5.353961 x 10~ 6 


2/d cos 2o\d 




-1.571576 x 10~ 6 
1.3568827 x 10~ 5 
1.7333956 x 10~ 4 . 



In Table 2, a — h are combinations of Laplace coefficients (e.g. Murray and Dcrmott 1999), and 



k E = H E f3 E m E , k D = HDP D m D,k 12 = Gm E m D HD0D'-> 



T 



j2c 



4.5Ci - 6C 2 , T j4c = 8C 3 - 10C 4 , T j2d = 3Ci - 4.5C 2 , T j4d = 5C 3 - 28C 4 ; 



Cll 



cuf3 E (J, E 

1 



c 2 



Cl 2 P D ^D 



r 7 

1 D 



c 3 



r ll 

1 E 



-GMm E R e J 2 ,ci2 



1 



-GMm D R e J2,C2i 



3 



C22/?£H 



pll 
1 D 



GMm E Rl.h,C22 



-GMm D RUi. 



2 " e " - 2 

2.2 5 Parameter. Tidal parameters and the rate of evolution through resonance 

The passage of the pair E-D through the 2:1 resonance is modeled here by including the time variation 
of the Hamiltonian's coefficients, following the same strategy given in (TW1988; TW1990). Due to tidal 
dissipative processes, the semi-major axes and eccentricities of the satellites suffer slow variations. Among 
all coefficients in Eq. ([2]), the most sensitive one to such variations is A. In fact, expanding the expression 
of A in a double series of up to first order, it can be shown that, in the proximity of resonance and for 
near-circular orbits, 



A » -(2n D - n E ) - AB{I E + I D ). 



(3) 
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The first member at right-hand side of Eq. ([3]) corresponds to the resonant combination of the mean 
motions which, in the proximity of the exact 2:1 commensurability, suffers large variations. All other 
coefficients of the Hamiltonian will be kept constant because they suffer small variations in the vicinity 
of the resonance. 

In order to study the system considering a wide variety of initial configurations in the vicinity of the 
exact resonance, we define the parameter 



S = 4A + 2(C + D). (4) 

In a first approximation and for small eccentricities, 6 ~ 4ri£) — 2n E — vj E — wd- The value of 8 changes 
sign in the middle of the exact resonance when we consider only the resonant term factored by K in 
Hamiltonian f2}. 

Adopting S as a free parameter, we can model the passage of our dynamical systems through the 
resonance. The rate of variation of S is given by 



8 « 4i (5) 
Axi±(2h D -h E )-4B(i E +i D ) (6) 



3 rii e 

In Eq. ([5]), the variations of the coefficients C and D have been neglected since their values (~ 
2 x 10 -16 ) are hundreds of times smaller than A. In order to evaluate A, we substitute the current values 
of the parameters of the system given in Table I in the following expressions: 



. , 63 n\k 2l M fRj\ 5 a 9k 2smi f R e V° 2 / 51 2 \ 

21 k 2l M fRA 5 57k 2smi (R e \ 5 

&i = — riiCi H — riiei, (9) 

2 Qi mi \ a,i ) 8 Qs M \a t J 

where, as usual, k 2 i, Qi, Ri refer to Love number, tidal dissipation function and mean radii of the 
satellites, respectively; the subscript S is used for Saturn (see Table I). Eqs. ([51 give us the time 
variation of the mean orbital elements of a satellite in a tidal model where the following approximations 
and hypothesis have been assumed (Ferraz-Mello et al. 2008): i) All dissipation functions are independent 
of the frequency of the corresponding tidal wave; ii) the satellite is in synchronous or in super-synchronous 
rotational state; iii) the satellite's mean motion is smaller than the frequency of the planet rotation; iv) 
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the tidal potential is truncated in second order of eccentricities, and the inclination terms have been 
neglected. The first and second terms on the right in Eqs. ([HI O are due to tides raised on the satellite 
by the planet, and tides raised on the planet by the satellite, respectively. 

The values of kn and Qi are unknown, and only some estimative of the domain of their plausible values 
are available. For instance, for Saturn, Qs ~ 18, 000 is the minimum possible value, which is obtained 
from the current configuration of Mimas (e.g. Murray and Dermott 1999). As the Enceladus-Dionc 
pair is currently in 2:1 resonance, we can roughly estimate the maximum value of Qs: if we consider 
that the pair E-D crossed the 2:1 resonance in the beginning of the solar system, then we find that the 
corresponding Qs must be less than 350,000; otherwise, the 2:1 resonance could not be reached during 
the age of the solar system. This estimative is based on the linear tide theory neglecting the terms of 
dissipation on the satellites in Eq. {SJ, and taking into account only E-D pair. Of course, a better 
estimative must be in close agreement with other resonant pair of the Saturn's system (Titan-Hyperion, 
Mimas- Tethys etc). For instance, assuming that the pair Mimas-Dione has crossed the 3:1 mean-motion 
resonance, Qs is given in the interval [18,000; 100,000] (Meyer and Wisdom 2008). 

In our numerical experiments we have used the intermediary value of Qs — 34,000. This value 
is slightly larger than that given by Ferraz-Mello and Hussmann (2005). Adopting Qs — 34,000, our 
numerical experiments show a very smooth evolution of the pair of satellites through the resonance, a 
condition which must be satisfied in order to avoid unrealistic evolutionary scenarios (TW1988). 

Fig. 1(a) shows A as functions of Qs. In our simulations, we have used A « 4.4x 10~ 14 , corresponding 
to the vertical full line in Fig. 1(a) (Qs = 34,000). The dashed line in Fig. 1(a) indicates the larger 
value for A, which is about 8.2 x 10~ 14 . 

2.3 Enceladus' Equilibrium eccentricity 

Consider a planet and one satellite where Eqs. (0 [9]) hold. Following TW1990, we calculate the non- 
resonant equilibrium eccentricity, for which the dissipation due to the tides on a synchronous or stationary 
rotating satellite ceases the growth of the semi-major axis due to planetary tides: 



~ / 1 n k 2l Q s {M\ 2 (RA 5 

e -"\/7^' D= k^sQ-{mJ UJ ' (10) 

Eq. (fl"U|) is obtained by setting n = in Eq. ©. In order to evolve through resonance, the satellite 
initial eccentricity must be less than the equilibrium eccentricity. The equilibrium eccentricity depends 
on the Love number and dissipation function of the satellites and the planet. In Fig. 1(b) we plot e eq . 
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of Enceladus as a function of k^E for some values of Qs, Qe- The two vertical lines at kiE ~ 0.0007 
and w 0.0018 refer to two values of the Love number for Enceladus found in literature (see Table I). No 
matter Qs, Qe and &2_e we adopt in Fig. 1(b), the important point here is that equilibrium eccentricity 
is e eq . > 0.02. Therefore we kept the initial value of the Enceladus eccentricity below the critical value 
0.02 in almost all numerical experiments shown in Section 3. 

2.4 Approaching the resonance zone 

The dynamics of the system @ depends strongly on the parameters (semi-major axes, masses, tidal 
constants), and on the initial conditions (eccentricities and initial critical angles). Following TW1988, 
we study the passage through the 2:1 resonance considering different initial values of eccentricities {fiEOi 
£do), keeping in all cases the same initial values of the semi-major axes (Table I). 

For each pair of initial eccentricities, we generate sets of initial conditions in the following way: without 
considering tide and adopting initially <je = &d = 0, Hamiltonian |J5J| is integrated over a short time 
span (about some years). Then we select on this orbit 20 different points, all of them having the same 
energy and slightly different eccentricities, but distinct values of ge and op. Each of these points will be 
used as the initial condition to be integrated for a long time, but now the tidal effect will be considered 
through the integration of 5. Each set with 20 numerical simulation will be denoted by run 1, 2, 3 etc in 
next section. 

3. Numerical Simulations 

In this section we describe the main regimes of motion of the pair E-D around the resonance considering 
different ranges of initial conditions, according to the methodology given in Section 2.4. After a detailed 
analysis of the main regimes of motion, we also estimate the probability of captures into some regimes 
of motion. For these tasks, we investigate the values of the eccentricities (e_E, e_o) and critical angles 
(<te, <td) reached by the orbits during the evolution through resonance. In order to better visualize the 
results, the plots of the variables are given in function of the S parameter. Let the total length of the 
time variation of S be divided in successive intervals A5 (typically A<5 « 10~ 5 or less). For each interval, 
in general we plot only the maximum and minimum value of the eccentricities and critical angles. 

Very far from the resonance encounter (5 -C 0), the long-term evolution of the pair of satellites is 
dominated by secular interactions, where the mean eccentricities oscillate regularly and the critical angles 
<je and on circulate (TW1988; Callegari et al. 2006). As the system approaches the resonance, different 
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scenarios of capture and escape from resonance may occur. 

3.1 Scenario for initial eccentricities smaller than 0.0035 

Let us begin considering values of initial eccentricities in two subintervals: < eEo ~ e D0 < 0.0015 
and 0.0015 < &eq ~ cdo < 0.0035. The main characteristics of the dynamics of the system inside the 
resonance in these intervals of initial eccentricity are summarized in Figs. 2(a-h). 

For < e^o ~ &do < 0.0015, erg is always captured into libration about zero in the beginning of the 
simulation. Fig. 2(c) shows an example where the capture occurs at S rs —0.015. Fig. 2(d) shows the 
evolution of <jd corresponding to the same orbit given in Fig. 2(c). <jd is captured slightly later into 
libration about 7r at S w —0.006. In this example, in spite of capture of <je, Enceladus eccentricity keeps 
its value near the initial one (Fig. 2(a)), as the system evolves through resonance. In the case of ctd, 
its capture implies a large growth of mean eccentricity of Dione (Fig. 2(b)). Based on previous studies 
(CY2007), the libration of oe resembles quite well the paradoxical (or small-eccentricity) libration. On 
the other hand, the libration of <td corresponds to a true libration (regime RDI in CY2007). We denote 
the aforementioned motions of <7e and ue by PL(0) (paradoxical) and L(7r) (true), respectively. 

For S > —0.006, Figs. 2(c,d) also show that the libration of <te, <jd and Aw occur simultaneously. 
This kind of motion of the system was first described in CY2007, and may occur when the system attains 
values of energy near the maximum. In order to better understand this regime of motion inside resonance 
(i.e, simultaneous libration of the critical angles), we compute surfaces of sections at several values of S 
(Fig. 3). Enceladus sections are defined through the condition yo = (see Eq. ([2])), represented on the 
plane {ceCOSGe X sin erg). Dione section is defined through the condition %/e = 0, and is represented 
on the plane (eccoscfl x e^smCD) (CY2007). Figs. 3(a,b) show the section for 5 = +0.007. The 
energy is H = 18.2 x 10~ 12 , very near the maximum of energy (H w 18.22 x 10 -12 ). In this portion of 
the phase space, the system is in the co-rotation zone, where libration of the critical angles associated to 
resonance and long-period oscillation of Aw can occur at same time. Fig. 3(c) shows the time variations 
of the elements of a single orbit where Aw and oe librate with same long-term period and amplitude. 
Fig. 3(d) shows the time evolution of the argument of pericenter of the orbits shown in Fig. 3(c): w\ are 
circulating with the same frequencies since Aw librates. 

In the interval of initial eccentricities 0.0015 < ceo ~ eoa < 0.0035, another type of motion often 
occurs: <je always circulates during the evolution through resonance (Fig. 2(g)), while ud is captured into 
libration about n (Fig. 2(h)). In this case the conjunctions between satellites occur in a line oscillating 
with large amplitude around the apocenter of Dione orbit. Note that Aw circulates (Fig. 2(h)). In this 
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example we can show that the system evolves in direction of maximum energy, but does not enter in the 
co-rotation zone (see Figure 11 in the Appendix). 

We can summarize our results on the evolution through resonance with initial eccentricities in the 
interval < eso ~ e D0 < 0.0035 as follows: in all 140 numerical simulations done considering this range 
of initial conditions, ctd is always captured into libration around n. In the subinterval e i0 < 0.0015, the 
dynamics is well characterized by the example given in Figs. 2(a-d). In the case where 0.0015 < eso < 
0.0035, the dominant regime of motion is given by the example shown in Figs. 2(e-h). 

We finish this section discussing the run 3, where initial eccentricities are of the same order of the 
current values (see Table I). For these eccentricities, the dominant regime of motion in almost all sim- 
ulations is similar to that given in Figs. 2(e-h). However, some exceptions may occur. For instance, 
Figs. 2(i-l) show a case where, after the system has entered in resonance without capture of oe (similar 
to Figs. 2(e-h)), a sudden reduction of eg occurs at S ~ (Fig. 2(i)), and <je is finally captured in 
the paradoxical libration PL(0) (Fig. 2(k)). At S +0.018, the amplitude of libration increases when a 
jump in eg appears. The same occurs with oe at <5 ~ +0.028, in this case, leading oe to return to the 
circulation regime. The jumps in the curve of eg occur since the system cross the separatrix associated 
to some secondary resonance (see also TW1990). 

Secondary resonances in the system E-D were studied in detail in CY2007, and they occur in portions 
of the phase space where some linear combination of oe and a d is commensurably. The centers of some 
secondary resonances in our system are illustrated in Figs. 3(e,f): the fixed points indicated by 1/1 and 
1/2 correspond to the loci in the phase space where oe + &d and 2a E + o D librate about it, respectively. 

Finally we note that the behavior of an in the example given in Figs. 2(i-l) is the same as that shown 
in Figs. 2(b,f), except that at S ~ +0.002 the amplitude of libration increases, while it seems to cause 
larger oscillation of er> (Fig. 2(j)). 

3.2 Scenario for 0.0035 < e E0 « e D0 < 0.017 

Figs. 4(a-h) show two examples taken from run 8, where e^o ~ 0.007. The system crosses the resonance 
and oe is captured into libration at S m +0.01. However, in opposition to the case shown in Figs. 2(a,c), 
the mean eccentricity of Enceladus now increases significantly during the evolution through resonance 
(Fig. 4(a)). This kind of motion corresponds to the capture of ge into true libration about zero, which 
will be referred as L(0) hereafter. The fixed point associated to L(0) and the chaotic separatrix associated 
to resonance can be seen in surface of section shown in Fig. 3(f). The center of the true libration is 
located at e^cosoB ~ +0.026 for S = +0.025. Note in Fig. 3(f) the presence of several high-order 
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secondary resonances near the separatrix of L(0). See also Figures 12(c,d) in the Appendix. 

In the case of Dione, its mean eccentricity (Fig. 4(b)) suddenly decreases at 5 w +0.003, showing 
that the system passes through resonance without capture into libration of gd- Fig. 4(d) shows gd 
circulating during the evolution through resonance. 

The case shown in Figs. 4(e-h) is similar to the example discussed in Figs. 4(a-d), except that now, 
for 8 w +0.016, the amplitude of libration of ge increases (Fig. 4(g)). Whenever ge increases for a given 
5, we have also a corresponding effect on e# and en (Figs. 4(c,f)). Large variation in the amplitude of 
the critical angle can cause the system escapes for larger values of S. See also Meyer and Wisdom (2008). 
As we will see in Section 3.3, escape from resonance by natural increasing of amplitude of critical angle 
is very often for larger values of initial eccentricities. 

3.3 Scenario for e E0 ^> e D0 

Consider now the results of numerical simulations where the initial conditions are such that 0.0035 < 
e_oo < 0.007 and eEo > 0.017. As we will see in Section 4, in these ranges of initial eccentricities the 
probability of capture into libration of <je around zero is small. Figs. 4(i,k) show an example where the 
system passes through resonance without capture of ge- In the trajectory, gd is temporarily captured 
into libration about 7r at S w (Fig. 4(1)) since the system escapes from resonance at 6 w +0.016 by 
increasing of the libration amplitude of gu (see discussion in Section 3.2). 

Figs. 5(a-d) show an example where the system crosses the resonance without any capture. In fact 
only a very short temporary libration around zero at 5 w +0.015 can be seen in Fig. 5(a), but the 
evolution of oe can be considered very similar to the case shown in Fig. 4(i). In the case of Dione critical 
angle, gjj is not captured (Fig. 5(b)). On the other hand, in the trajectory shown in Figs. 5(e-h), the 
escape of the system from L(0) regime of motion occurs after a long period of temporary capture. 

Here, let us discuss some complex evolutionary scenarios of capture and escape from secondary reso- 
nances which appear in several intervals of initial eccentricities discussed in this work. This time, none 
of ue or ur> is captured into their main libration centers (Figs. 6(c,d)) but, in opposition to the cases 
shown in Figs. 5(a-d), both eccentricities increase (Figs. 6(a,b)). The growth in eccentricities occurs due 
to some capture into a secondary resonance (see Section 3.1). Figs. 6(a-d) show examples where capture 
into ge + o\d secondary resonance occurs at 6 w +0.005: the angle ge+gd librates while ge and gd 
circulate. The scattered points at 5 w +0.025 in Figs. 6(a-d) belong to a chaotic region of the phase 
space which can also be identified in some surfaces of section, more precisely, in the Enceladus section 
(Fig. 3(f)): the source of this chaos is the separatrix associated to the true libration of ge (regime L(0)). 
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Other examples of the growth of Enceladus eccentricity due to capture into secondary resonances are 
given in Figs. 6(e-l). This time, we show the libration of oe + &d in Fig. 6(h) and <je + 4cr£> in Fig. 
6(1). In both cases, the system escapes from secondary resonances at large values of 8. This occurs due 
to interaction of their domain with the separatrix of true resonance L(0). Fig. 6(g) shows that, after 
the escape, oe is briefly captured into libration and escapes. In the case shown in Fig. 6(k) we see an 
interesting feature: after the escape from secondary resonance, oe is captured into L(0) with decreasing 
amplitude while e# increases. In the case of Figs. 6(i-l) it is worth to note how high order secondary 
resonance is able to excite the initial eccentricities of the satellites. 

The role of secondary resonances in the dynamical history of natural satellites was studied by several 
authors. As far as we are aware, capture into secondary resonances without prior libration of the system 
into libration of <je or op, as shown in Fig. 6, seems to be new, not described in previous works (TW1989; 
TW1990; Malhotra and Dcrmott 1990; Henrard and Sato 1990). 

To close this section, it is important to mention some differences when comparing the general three- 
body model with simplified models of resonance: in Fig. 4(i) and Figs. 5(a,b,f), the final mean eccentricity 
after passage through the resonance is smaller than its initial value. This is the general behavior for one- 
degree-of-freedom models (see Peale 1986), in the case of passage without capture. However in Fig. 4(j), 
Figs. 6(e,f,j) and Fig. 5(e), the final mean eccentricity of the satellite is larger than the initial one, which 
is not predicted by integrable models (TW1988). 

3.4 Scenario for e E o <C e D0 

Since the Dione mass is about 11.2 times larger than the Enceladus one (Table I), Dione's perturbation on 
Enceladus is stronger than the Enceladus' one over Dione. Thus, the hypothesis to consider e^o <C e^a 
may not be completely realistic unless some previous phenomenon outside the 2:1 resonance might have 
happened, pumping en to higher values. In this section we only show the results for ceo <C cdo, leaving 
a brief discussion on the reliability of this initial conditions to Section 4. 

We already showed that quite different regimes of motion are found when we consider e#o ~ e^o 
(Sections 3.1 and 3.2), and e# S> (Section 3.3). For e^o <C e^o, new regimes of motion arc 
obtained. Figs. 7(a-d) show a case where eoo ~ 0.004 and e#o « 0.002. At 6 w —0.002, <je is temporary 
captured into (Fig. 7(c)), and librates until 5 w +0.004, when the system escapes from resonance. After 
the escape, however, Enceladus and Dione eccentricities increase suddenly, since the system is captured 
into a secondary resonance, where the angle 2a E + ®d librates around 7r (Fig. 7(d)). At S w +0.036 the 
system escapes from secondary resonance. 
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The mechanism of capture into secondary resonance shown in Figs. 7(a-d) is different from those 
discussed in Fig. 6: there, the system was automatically captured into secondary resonances, while 
now the system is first captured into the main libration, escapes, and enters in a secondary resonance. 
This mechanism of capture into secondary resonances is similar to those reported in TW1989, TW1990, 
Malhotra and Dermott (1990). 

Figs. 7(e-h) differ from the previous case (Fig. 7(a-d)) since <je is sooner captured into libration about 
0. At S w —0.004, the system escapes from libration (Fig. 7(g)), and oe circulates until 5 w +0.005, 
when a new temporary capture into libration about zero occurs. The system escapes and enters into the 
(7 e + <j d resonance (Fig. 7(h)). 

We finish Section 3 showing two examples of complex regimes of motion found only when ceo <C eoo 
and initial eccentricities are very near zero. In the case shown in Fig. 7(j), <7e is captured in paradoxical 
libration in the beginning of the simulation, but escapes and is again re-captured at S w —0.005. Note 
the similarity between Figs. 2(a) and 7(i). The case shown in Fig. 7(1) is more complicated since the 
system is alternating between libration and circulation during the evolution through the 2:1 resonance. 

3.5 Statistics 

In this work a total of at least one thousand numerical simulations were performed. Although not 
exhaustive, we think that a significant part of the resonant domain could be covered. These simulations 
revealed several different dynamics, and most of them can be collected in one of the following cases: 

1) Capture of ue into paradoxical libration about zero (PL(0)), and capture of ctd about tt (L(tt)) 
(Figs. 2(a-d)). Sometimes, capture of oe occurs at large values of 6 and, eventually, escapes (Figs. 2(i-l)). 

2) The same as 1) where <7e alternates between libration in paradoxical libration and circulation (Figs. 
7(i-l)). 

3) Passage through resonance without capture of o\e, and capture of an into libration about tt (Figs. 
2(e-f)). 

4) Capture of <je into true libration about zero (L(0)), and passage through resonance without capture 
of an (Fig. 4(a-h)). 

5) Escape from temporary libration of <je in L(0), and passage through resonance without capture of 
a D (Figs. 5(e-h)). 

6) Passage through resonance without capture of oe and an- Sometimes oe quickly librates about 
zero (Figs. 5(a-d)). 
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7) Passage through resonance without capture of erg into L(0), and capture of ctd into L(ir) with 
subsequent escape (Figs. 4(i-l)). 

8) Automatic capture into secondary resonances (Fig. 6). 

9) Temporary capture of ge into L(0), followed by capture into secondary resonance (Figs. 7(a-d)). 
In some cases, the system is re-captured into L(0) at large values of S. 

10) The same as 9, except that ge is captured into libration in the beginning of the simulation (Figs. 
7(e-h)). 

11) Exceptions corresponding to cases which do not match any case listed above. 

Fig. 8(a) shows that, for initial eccentricities < do < 0.0015, only the case 1) is possible. For 
0.0015 < eto < 0.0035, ge is not captured into paradoxical libration (except the three cases when 
e^o = 0.0025), while the capture of gd about 7r is certain and the most probable event is the case 3). 

For initial eccentricities taken in the interval 0.0035 < e^o < 0.016, the capture into state 4 is a very 
probable event. By inspection of Fig. 8(a) we see that the probability capture of g e into libration about 
zero increases (slightly) as long as the initial eccentricities are selected in the interval 0.005 < e^o < 0.01. 
This result is in opposition to those predicted by single resonance models, which states that capture 
into true libration always occurs for initial eccentricities smaller than some value e ma x (see discussion in 
Section 4). The cases 6 and 7 never occur for 0.0035 < e i0 < 0.016, and the case 5 is a very rare event 
in this range of initial eccentricities. However, there are some cases of escape from secondary resonances 
(e.g. Fig. 6(c)), resulting in a slight increase in the total number of escapes for 0.0035 < e i0 < 0.016. 

Fig. 8(b) shows the statistics for e E o ~> e Da . In the two runs where e^o = 0.008, e D0 = 0.004, 
and eEo = 0.01, e^o = 0.005, the dominant regime of motion is L(0). In all cases where eEo > 0.017 
and e_oo > 0.0035, the probability of passage through the resonance without any capture is high. When 
e_Eo > 0.017 and eoo < 0.0035 (as the case eEo = 0.02, euo = 0.002), although capture of ge in general is 
ruled out, gu always falls into libration in 7r. Other interesting run is the case e^o = 0.0045, eoo = 0.0022 
corresponding to the values of the current eccentricity: the dominant regime of motion is the case 3, and 
there are two exceptions where paradoxical libration appears for large values of 5 (Figs. 2(i-l)). 

Fig. 8(c) corresponds to the cases e^o ^ e Do- I n this range of initial eccentricity, the new features 
are the presence of: i) capture into secondary resonances according to the cases 9 and 10; ii) several 
complex scenarios of capture and escape when the initial eccentricities are both very near zero (note the 
cases eEo = 0, eoo = 0.002 and e_Eo = 0, eoo = 0.004). 
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3.6 Parametric Plane 

From Eq. @ we write A = | [6 — 2(C + D)], so that the Hamiltonian © becomes a function of S. 
H and S are the natural parameters to describe the evolution of the system around the resonance, as 
pointed out in Section 2.4. In this section we show the parametric plane (<5, H), which is very convenient 
to summarize the results about the dynamics shown in previous subsections. 

The parametric plane (Fig. 9(a)) is constructed by calculating the zero-gradient points of Eq. ([21), in 
the case of symmetric configurations (i.e., ue, &d = 0, 7r). Zero-gradient points can be seen in Fig. 9(b) 
in the representative plane of initial conditions, which are given by the level curves of the H = H(eE, &d) 
for a fixed value of 6. In Fig. 9(b) shows also, in bold lines and gray symbols, the curves defined by &e = 
and <td = 0, respectively, which must intercept themselves in the zero-gradient points (e.g. CY2007; see 
also Figure 10 in the Appendix). 

For (j e , <Jd = 0,7r, Eq. ([2|) has at most five zero-gradient points: C, D, M, O and P. In the plane 
(S, H), the energy of each of these points as function of 5 generate five curves. Fig. 9(a) shows the curves, 
which have been labeled with the same letter of the corresponding critical point. 

The main characteristics of the resonant dynamics can be drawn from the curves plotted on the 
parametric plane. Very far from the resonance and for negative values of (5, only the critical point M 
exists. In this case, the dynamics of the system is dominated by secular interactions (TW1988). 

For 5 as —0.005 a bifurcation occurs in the parametric plane. For —0.005 < 6 < +0.0088, only the 
curves C, D, M delimit different regions of initial conditions. The current position of the pair E-D in the 
parametric plane belongs to the region between C and D, and the corresponding dynamics of the system 
was studied in detail by CY2007. See also Figure 12 in the Appendix. Between D and M we have the 
co-rotation zone (see Section 3.1). For energies larger then M, no motion is possible, and we have the 
forbidden regions of initial parameters. 

For 6 > +0.0088, a new bifurcation occurs in the plane (S, H), and we can see the rise of the points O 
and P. There are important regimes of motion in this range of (5, as described in Section 3.2. For instance, 
the rise of the separatrix of the true libration of oe around zero (L(0)), may occur for S > +0.0088 (see 
Fig. 3(f)). 

The evolution of the trajectories shown in Section 3 can also be shown in parametric plane (see 
TW1988). In fact, as tidal evolution changes the values of S and H, the system is driven to different 
regions of the parametric space. In Figure 11 in the Appendix we show some examples of trajectories 
evolving in the plane (S,H). 
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4. Discussion 



In this section, we review the scenarios of the evolution of the pair E-D through 2:1 resonance given by 
intcgrable models of resonance, and compare them with our results (Section 3). The formation of the 
current configuration of the pair E-D is also discussed. 

In a scenario of tidal evolution, Sinclair (1972) and Henrard and Lemaitre (1983) found that the prob- 
ability that the pair E-D has passed through 2:1 resonance without capture of <jd around ir is about 81% 
for eoo > 0.0032. They considered e# = and neglected long-period terms, that is, they used simplified 
models of the so-called "Dione resonance" . Remember that in adiabatic evolution through resonances 
considering one-degree-of-frecdom models, capture into libration is certain only for initial eccentricities 
below a determined value e max . Sinclair (1972) proposes that, before the resonance encounter, Dione 
initial eccentricity was of the order of e max ~ 0.0032, slightly larger than the current value (ssO.0022), 
and <jd was circulating in retrograde direction. The passage through resonance without capture reduced 
e_o to values near the current one, and changed the direction of circulation of an- 

According to our results, capture of <jd into true libration is a certain phenomenon whenever eno < 
0.0035 (see Sections 3.1 and Fig. 8). Moreover, the probability of the system to pass through resonance 
without capture of <r D is high for e Da > 0.0035. We conclude therefore that, in the case of the "Dione 
resonance" , our results are in good agreement with the classical ones. 

In the case of integrable models of the "Enceladus resonance" , where eo = and all long-period 
are neglected, classical results are the following: capture of oe into libration about zero is guaranteed 
provided that e_eo < e max - Peale (1986) and Sinclair (1972) give e max ~ 0.017 and w 0.019, respectively. 
Considering the current data given in Table I, and using the equation of e max given by Peale (1986), we 
obtain e max w 0.0183. As currently e# w 0.0045, i.e., e# <C e max , classical results consider that capture in 
libration was automatic in the past. The current small values of the free Enceladus eccentricity (<~ 10~ 4 ; 
see Ferraz-Mello 1985), and amplitude of libration (1.505°), would be related to tidal dissipation within 
the satellite. 

The results shown in Section 3.2, where 0.0035 < do < 0.017, show that the probability of capture 
of g e into true libration (L(0)) is high, in good agreement with classical results. However, in the case of 
capture into L(0), the Enceladus eccentricity could have increased to values larger than the current one. 
Since natural escape from this libration (case 5 in Section 3.5) is a very rare event (see Fig. 8(a)), the 
system could not have evolved to the current configuration. 

In Section 3.3 we have studied the dynamics for c_eo ^> cdo- In the cases where e^o is somewhat 



16 



(> 0.017) larger than the current one, the number of captures into L(0) is small, and the number of 
unsuccessful passages and escapes are large, in good agreement with the results obtained with integrable 
models resonance. Numerical experiments with eEo > 0.017 and 0.0035 < eoo < 0.007 have been done in 
order to find examples where and <je could be reduced to the current values when the system crossed 
the 2:1 resonance. In spite of several tests, none of the experiments seems to be able to reproduce the 
current configuration of the E-D system. 

Other scenarios obtained from our simulations, not discussed in classical works, are related to capture 
of the system into paradoxical (small-eccentricity) libration of <je about zero. In Section 3.1 we have 
shown that, in order to get capture into paradoxical libration, we must have ao < 0.0015. As we deal with 
general three-body problem taking into account the mutual interaction between E-D, for e^o < 0.0015, 
<jd usually falls into L(7r) regime of motion (see discussion above). This is in complete opposition to what 
we would expect considering the current configuration of E-D system. The capture into L(7r) excites the 
initial eccentricity of Dione, so that this libration state of ctd and resulting eo are not compatible with 
the current scenario of the E-D system. 

The last case considered in Section 3 corresponds to the case where the Dione eccentricity is larger 
than the Enceladus one (Section 3.4). According to our results discussed above, if eno > 0.0035, capture 
of (To into 7r is avoided and, if e^o < 0.0015, capture into small-eccentricity libration is certain. These 
two simultaneous events would be enough to recover the current configuration of the system. Some cases 
corresponding to the runs where e E — (e.g. Figs. 7(e-h)) agree with this hypothesis: in the beginning 
of the integration, a E is always captured into paradoxical libration while ctd circulates. However, the 
subsequent evolution of the system is marked by a complex motion, similar to the case 10 discussed in 
Section 3.5, where the system enters into a web of secondary resonance. 

At first glance, the hypothesis eno eEo is not plausible. Meyer and Wisdom (2008) also show some 
simulations of the evolution through resonance of the pair E-D where eoo 3> &eo- They also propose 
some additional mechanism to increase the Dione eccentricity, for instance, the passage of Dione and 
Mimas through the 3:1 mean-motion resonance. Alternatively, although not so probable, we can also 
think that some increase of the eccentricities of Dione or Enceladus could have happened during the 
planetary migration (Tsiganis et al. 2005; Nogueira 2008), if giant planets as Uranus and Neptune had 
approached very close to Saturn. 

In our model, we do not reproduce the current configuration of the pair E-D in the case where 
euo 3> e Eo since Enceladus eccentricity in general increases, due to capture of the system into secondary 
resonances, reaching values much larger than current 0.0045. In a similar scenario, Meyer and Wisdom 
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(2008) show an example where eccentricities decay and they can reproduce, in a limited time interval, 
the current configuration of the pair E-D. However, there is a main difference between the methodology 
applied by Meyer and Wisdom (2008) and ours. As pointed by TW1988, tidal friction acts in two different 
forms: the time variation of the Hamiltonian's coefficients, and the direct effects of tides on the degrees 
of the freedom of the problem. The latter was considered by Meyer and Wisdom (2008), but not in this 
work. Maybe an extension of our model by including satellites and planetary tides directly in canonical 
equations could lead to more exciting scenarios of evolution of the Saturn-E-D system. 

5. Conclusions 

We have studied the dynamics of the pair E-D during the passage through the 2:1 mean motion resonance 
with a two-degrees-of-freedom planar model. We have collected the essential physics of the general planar 
three-body problem (TW1988, CY2007), so that our basic dynamics is a generalization of some simplified 
models used by Sinclair (1972), Henrard and Lemaitrc (1983) etc. 

Under the action of tides, the system is slowly evolving toward the resonance, starting from several 
different initial conditions. The choice of these initial conditions is mostly dictated by the possible 
eccentricities of both satellites in the past. The simplified theory of the equilibrium eccentricity for a 
synchronous satellite is used to obtain a rough estimate of the range of the eccentricities. 

Based on this estimative of initial eccentricities, hundreds of initial conditions (almost 1000) have 
been integrated considering slow tidal evolution with Qs — 34, 000, fes = 0.34. We summarize our main 
results, and the main regimes of motion associated to 2:1 mean-motion resonance are: 

•small-eccentricity (paradoxical) libration and true libration of Enceladus critical angle (<Je), both 
about zero; 

•true libration of Dione critical angle (<7d) about 7r; 

•libration into several secondary resonances involving linear combinations of oe and oo- 
All these regimes of motion appear in sets of runs where different initial conditions have been consid- 
ered: 

•For initial eccentricities smaller than 0.0035, ur> is always captured into libration about n, so that 
Dionc eccentricity (eu) increases to values larger than the current one. Moreover, in order to guarantee 
the capture of <je into small-eccentricity libration (similar to the current one), our model requires ceo < 
0.0015. The simultaneous occurrence of both events conflicts with the current scenario of E-D system. 

•For 0.0035 < e i0 < 0.017, we have found that there is high probability of capture of oe into true 
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libration about zero and, as a consequence, be increases to values much larger than the current one. 

•In the case where e^o > 0.017 and e^o > 0.0035, the number of captures of erg, &d into libration is 
small. However, in several runs, capture into secondary resonances increases eg. 

•When eoo S> ceo, the determinant role in the evolution of the system is played by secondary 
resonances and complex regimes of motion involving exchange of different types of captures. In spite of 
the rich dynamics, no example evolving to the current E-D configuration has been verified. 

We conclude that the current libration in E-D system cannot be reproduced with over-simplified 
models of resonance. Even with our two-degrees-of-freedom model, where several complex and chaotic 
motion are possible, we cannot recover in high details the current configuration of the system. More 
general models including more satellites and tidal dissipation in all bodies, seem to be necessary to 
recover the Enceladus-Dione history. 
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Appendix 
Figure 10 

In Fig. 9(b)) we have shown levels curves of H — H{eE,e-D) for a fixed value of 5. Fig. 10 shows 
additional examples of representative planes of initial conditions for different values of 5. In Figure 10 
we have also plotted, in gray and bold lines, respectively, the curves defined by &e — and &jj = 0. 
As discussed in Section 3.6, the evolution of the level curves can give important information about the 
dynamics of the system. 

Figure 11 

As pointed in Section 3, as the system evolves due to the tidal dissipation, the values of 5 and H change, 
and the system passes through different regions of the parametric plane defined in Fig. 9(a). In Fig. 11, 
we show the variation of the energy of some trajectories in the parametric plane. Trajectories starting with 
eccentricities eio < 0.0015 (e.g. dashed curves in Fig. 11(a)), in general enter in the co-rotation zone, the 
region located region between the curves M and D where both critical angles can librate simultaneously 
(see Section 3.1). For slightly larger values of eccentricity, the system evolves to the regions between D 
and C, far from co-rotation zone. One example is shown in Fig. ll(a,b) by gray line, and corresponds to 
a trajectory of case 3) discussed in Section 3.5. For initial eccentricities e i0 > 0.0035, the system is, in 
general, captured into L(0) or into some secondary resonances. The evolution of such type of trajectories 
are illustrated in Fig. 11(b) by dashed curves. 
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Figure 12 

In order to help the understanding of the main regimes of motion attained by the system during its 
evolution on the parametric plane, we construct several dynamical maps at frozen values of S (Fig. 12). 
They are based on the spectral number N (see Ferraz-Mello et al. 2005 and references therein), and 
are built integrating the averaged equations of motion in a grid of points in the representative plane for 
a given S. For each point, we obtain numerically the Fourier spectrum of the solution. The spectral 
number N is defined as the number of significant spectral peaks which amplitude is greater than 0.1% of 
the highest peak of the variable x E w e# cos erg. 

In Fig. 12, we see the periodic orbits associated to the main regimes of motion of the critical angles 
indicated in Fig. 9(a): paradoxical libration of oe around zero (denoted by PL(0)); true libration of 
oe around zero (L(0)); true libration of <jd around n (L(tt)); paradoxical libration of ctd around zero 
(RDIP). In order to better visualize the loci of PL(0) and L(0), we also plot in Figure 12 the loci of 
&e = (see also Fig. 9(b)). The black regions seen in Fig. 12 correspond to the separatrix of the regime 
L(7r) associated to Dionc resonance. The oval structures are associated to secondary resonances between 
<je, Near the maximum value of energy (M) we have the co-rotation zone. 

Figures 12(a,b) are similar to Fig. 8 given in CY2007, where the dynamics of the system was studied 
in details in these regions of the parametric plane. For larger values of S (Figs. 12(c,d)), however, the 
distribution of the periodic orbits and the dynamics of the system are quite different. The maximum 
forced Enceladus eccentricity associated to the periodic orbit PL(0) and Dione eccentricity associated 
to L(7r), increase in the parametric plane. For larger values of S, a new chaotic zone rises between the 
domain of secondary resonances, which is associated to the regime of true libration of oe around zero. 
In Figs. 12(c,d), we denote the periodic orbit associated to this regime of motion by L(0). Inspection of 
Fig. 12 shows that the loci of &e = now represents the center of true libration L(0). See also Fig. 3(f). 
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Captions 

FIGURE 1. (a) A as function of Qs- Vertical dashed line indicates an estimative of the minimum 
value of Qs- Vertical full line indicates the value of Qs used in the calculation of the 5. (b) Equilibrium 
eccentricity of Enceladus (e eq .) as function of k 2 E for several values of Qs and Qe- k 2 s = 0.34 (see 
Table I). From top to bottom: dashed curve: Qs — 18, 000, Qe — 100; dashed-bold curve: Qs = 34, 000, 
Qe = 100; full line: Q s = 18,000, Q E = 20; full-bold line: Q s = 34,000, Q E = 20. The vertical lines 
indicate two values of k 2 E found for Enceladus (see Table I). 

FIGURE 2. Maximum and minimum values of the eccentricities and critical angles in fixed intervals of 
AS, for different initial conditions, (a-d) e^o ~ e D0 ~ 0.001. (e-h) eso ~ e D0 ~ 0.002. (i-1) e#o ~ 0.0045, 
e_oo ~ 0.0022. In (d,h,l), the gray points show the evolution of Aw = ge — &d = — we- 

FIGURE 3. (a) Enceladus section for S = +0.007 and energy in the co-rotation zone around the 
maximum value M (see Figs. 12 and 13 in the Appendix), (b) Dione section corresponding to (a), (c) 
Time evolution of ct_e, gd and Aw for initial conditions in the co-rotation zone: e^o = 0.0016, geo = 0, 
Abo = 0, e/5o = 0.007856, gdo = n, Xdo = v. (d) Time evolution of the pericenter of Enceladus (full 
circles) and Dione (open circles) of the orbit shown in (c). (e) Enceladus section for S = +0.007 and 
energy far from co-rotation zone. The domain of some secondary resonances are indicated by 1/2 and 
1/1. (f) Enceladus section for S = +0.025. 

FIGURE 4. Maximum and minimum values of the eccentricities and critical angles in fixed intervals of 
AS, for different initial conditions, (a-h) Two distinct trajectories from run 8: e# ~ ~ 0.007. (a-d) 
geo = 35.573°, G Da = 153.475°. (c-h) geo = 1.456°, g do = 130.431°. (i-1) Evolution of one trajectory 
from run 7: e E a « 0.017, e D0 ~ 0.0035. 

FIGURE 5. Maximum and minimum values of the eccentricities and critical angles in fixed intervals 
of AS, for two distinct ((a-d) and (e-h)) trajectories from run 10, where e^o ~ 0.02, eno ~ 0.007. 

FIGURE 6. Maximum and minimum values of the eccentricities and critical angles in fixed intervals 
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of AS, for different initial conditions, (a-d) e_eo ~ 0.017, eoo ~ 0.0035. (e-1) e# ~ 0.02, eoo ~ 0.0035. 
Figs. 6(h,l) show the evolution of the angles <je + &d, o\e + 4<td, respectively. 

FIGURE 7. Maximum and minimum values of the eccentricities and critical angles in fixed intervals 
of AS, for trajectories with e^o "C £dq- ( a "d) e^o ~ 0.002, eoo ~ 0.004. (e-h) e_Eo ~ 0, e^o ~ 0.008. 
(i,j) The same as (a,c), respectively, where eEo ~ 0, e^o ~ 0.002. (k,l) Same as (a,c), respectively, where 
e E0 ss 0, epo ~ 0.004. Figs. 7(d,h) show the evolution of the angles 2tr B + <r D , oe + &d-> respectively. 

FIGURE 8(a,b). Number of occurrences of the regimes of motion 1-11) listed in Section 3.5 as function 
of the initial eccentricities. The final values of S attained in different runs are shown on the top of the 
bars, (a) e^o ~ er>o- Sfi na i ~ 0.0324 in the interval e^o = and e^o = 0.003. (b) e^o 3* e D0- (c) 

FIGURE 8(c). Continued. 

FIGURE 9. (a) The paths of zero-gradient points M, D, C, P and O in the parametric plane (5,H). 
The localization of the main regimes of motion are indicated by symbols; black triangle: retrograde or 
direct circulation of <je, od; big circle: paradoxical libration of <je, &d around zero; •: true libration of 
a E around zero; black square: true libration of ao around 7r; square: paradoxical libration of <je around 
7r; secondary resonances, (b) Representative plane (eE,eo) of initial conditions for 5 = +0.018. Each 
curve is a level curve of the Hamiltonian for symmetric configurations. Positive and negative values 
of eccentricity correspond to Oi = and n, respectively. The gray points and bold line gives the loci of 
&e = and &d = 0, respectively. 

FIGURE 10. Representative planes (eE,en) of initial conditions for different values of S parameter. 
Each curve is a level curve of the Hamiltonian (2), for symmetric configurations. Positive and negative 
values of eccentricity correspond to <7j = and tt, respectively. The gray and black-bold lines are the loci 
of & E = and & D = 0, respectively, (a) S = -0.014. (b) S = +0.005. (c) 5 = +0.0088. (d) S = +0.018. 
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FIGURE 11. The paths of zero-gradient points M, D, C, P and in the (S,H) plane, (a) Typical 
trajectories taken from runs with eccentricities near zero: Run 1 (e^o ~ e D0 ~ 0.001, dashed lines), Run 
2 (g£ ~ e Do ~ 0.002, gray line), (b) Several trajectories from Run 8 (e#o ~ e D0 ~ 0.007, dashed lines), 
and Run 2. 

FIGURE 12. Dynamical maps showing the spectral number TV, as defined in the text, for several 
values of 5. (a) S = 0. (b) 6 = +0.007. (c) 5 = +0.01. (d) S = +0.018. The bold-gray lines show the loci 
of &e = obtained analytically and match the loci of the periodic orbits of PL(0) and L(0). 
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